Nx = 256;

lc = 16;
dx = 1 / Nx;
epsilon = 1/Nx;

lambda = 1d-2;
lambdarnd = 1d-2;
lambdamod = 3d-1;



x = (0:Nx-1)'/Nx;

% Purely random
if(1)
  V = lambdarnd * rand( Nx, 1 );
end

% Filtered random noise
if(0)
  V = lambdarnd * rand( Nx, 1 );
  fV = fft(V);
  kcutoff = 30;
  ind = find(abs(k)>kcutoff);
  fV(ind)=0;
  V = real(ifft(fV));
  V = V / max(abs(V)) * lambda;
%  V = V * Nx / kcutoff;
end

% Periodic 
if(0)
  V = zeros(Nx,1);
  % V = lambda * cos( 2*pi / lc * (x/epsilon) );
  % mu = 1.904722371438797e-02;
  mu = 1.921471959676978e-02;
end

% Periodic + random
if(0)
  V = lambda * cos( 2*pi / lc * (x/epsilon) ) + lambdarnd * rand( Nx, 1 );
end

% Periodic + simple long wave mode
if (0)
  V = lambda * cos( 2*pi / lc * (x/epsilon) )  + ...
    lambdarnd * cos(2*pi * x );
end

% modulated periodic
if(0)
  V = lambda * cos( 2*pi / lc * (x/epsilon) ) .* (1 + ...
    lambdamod * cos(2*pi * x ) + lambdamod * sin(2*pi *2*x)); % +  0.01 ...
%      * randn(size(x)));
end

e = ones(Nx,1);
H = spdiags(0.5*[-e 2*e -e], -1:1, Nx, Nx) / dx^2  * epsilon^2; % Discretization of the laplacian
H(1,Nx) = -0.5/dx^2*epsilon^2;
H(Nx,1) = -0.5/dx^2*epsilon^2;

H = H + spdiags(V,0,Nx, Nx);

[EF, D] = eig(full(H));
D = diag(D);


